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Abstract 

|n the interest of conservation, the importance of having a large habitat available for a species is widely known. Here, we introduce a 
MH lattice-based model for a population and look at the importance of fluctuations as well as that of the population density, particularly 
Q with respect to Allee effects. We examine the model analytically and by Monte Carlo simulations and find that, while the size of 
■ ^ the habitat is important, there exists a critical population density below which extinction is assured. This has large consequences 
1^ with respect to conservation, especially in the design of habitats and for populations whose density has become small. In particular, 
^^e find that the probability of survival for small populations can be increased by a reduction in the size of the habitat and show 
' 'that there exists an optimal size reduction. 
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(N 1. Introduction 

. Extinction is becoming a greater and greater issue all over 
'"q the world and is a cause of extreme concern. It has been es- 

• ^ timated that anthropogenic extinctions are resulting in the 

loss of a few percent of the current world's biosphere, which 
g rts of magnitude 3 to 4 times the natural background rate 
jMavet al.l . [l995h . The World Conservation Union (lUCN), 

• ^ through its Species Survival Commission (SCC) develops 

Criteria to assess the extinction rate for plants and animals 
^ all over the world which enables them to keep a so-called 
■ ■ ■ Sed -List (www.redlist.org) of species which are threatened 
with extinction in order to promote their conservation. The 
list currently shows over 16,000 threatened species around 
the world - a 45% increase on the figure from the year 2000 . 

It has been shown by examining both discrete fescudero . 
[2005) and continuous (.Skel lam, 1951) populations that, at 
least analytically speaking, there exists a critical habitat 
size Lc above which survival of a population is assured. 
Here we examine what role the population density plays 
since, intuitively, one would expect that even for L > Lc, a, 
sufficiently large population would be needed for growth. 
Lattice bas ed mode ls are widely used in ecolo gy (see 



for example Tainakal . Il988t iDurrett and Levinl llQQSt 



Itoh et al.l . 120041 ) and so we introduce such a model that 
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incorporates birth, death and diffusion. Unli ke other sim- 
ilar models such as the contact process (e.g. Harrisl 1974 : 
"iool, here two individuals must meet in 



Obornv et al 



order to reproduce whereas one individual can die by itself. 
This results in negative growth rates for populations that 
fall below a critical density due to reproduction oppor- 
tunities becoming rare. In real populations, the positive 
correlation between size and per capita g rowth rate o f 
a population is known as the Allee effect (jAUeel . 1931 ). 
which has recently r eceived much interes t (e.g. to ennid . 
200i iHurford et atl . l2006t Ijohnson et ad] . |2006). If the 
Allee effect is strong enough, the population size may even 
decrease for small population sizes as in our model. Due 
to this behaviour, the effect has bee n examined w i th re- 
spect to extinction (see for example Amarasekar^ 1998 : 
Courchamp et al.1 . 19991 : [Stephens et al. . 19991 and ref- 



erences therein) but primarily deterministically and so 
without fluctuations in the population density. Since fluc- 
tuations are likely to be highly significant for small popu- 
lations, we include the effects of these by examining Monte 
Carlo (MC) simulations in the hope to gain a more realis- 
tic picture of the importance of the population density on 
the chances of survival. Further to the stochastic methods 
used to study Allee effects, such as stochastic differential 
equa tions (e.g.lDennisl,ll989l), discrete-time Markov-chains 
e.g. Allen et a l.. 2005) or diffusion processes (e.g. Dennid . 



20021) . lattice based models have space, as well as time, as 



a variable and take into account individuals, rather than 
just the macroscopic view of the population. 



After introducing the model in the next section, we exam- 
ine the AUee effects present in our model in Section 3, par- 
ticularly with respect to a sudden decrease in population 
in Section 5. The effects of the fluctuations are examined 
in Section 4. 

2. The Model 

We have a d-dimensional square lattice of linear length 
L where each site is either occupied by a single particle 
(1) or is empty (0). A site is chosen at random. If the site 
is occupied, the particle is removed with probability pd, 
leaving the site empty. If the particle does not die, a nearest 
neighbour site is randomly chosen. If the neighbouring site 
is empty, the particle moves to that site. If however the 
neighbouring site is occupied, with probability the 
particle reproduces, producing a new particle on another 
randomly selected neighbouring site, conditional on that 
chosen square being empty. We therefore have the following 
reactions for a particle A: 



Act) < — > 't'A A + A 



3A and A 



(2.1) 



where represents an empty site. A time step is defined 
as the number of lattice sites and so is equal to approx- 
imately one update per site. We use nearest neighbours 
and, throughout most of the paper, periodic boundary con- 
ditions which, although more unrealistic than, say, reflec- 
tive boundary conditions, allow for better comparison with 
analytical results, since periodic systems remain homoge- 
nous. We later, however, examine some results with reflec- 
tive boundary conditions. 

Due to the conflict between the growth and decay pro- 
cesses in the model, we expect that with certain values of 
Pb and pd, extinction of the population would occur. In- 
dee d, many models d i splay i ng such a conflict (see for exam- 



ple 'VesDignani et al. V200(^;' Dammer and Hin richseiil . l2003l : 
lObo rny et al.. . ,2005. ; Peters and Ncclin, 2006) show a criti- 
cal parameter value separating an active state and an inac- 
tive or absorbing state that, once reached, the system can- 
not leave. As the rate of decay increases, the so-called order 
parameter (often the density of active sites) decreases, be- 
coming zero at a critical point, marking a change in phase 
or phase transition. In our case, the absorbing state would 
represent an empty lattice and so extinction of the popu- 
lation. 

To show that this is indeed the case for our model, we de- 
rive a so-called mean field equation (e.g. lOpper and Saadl . 



20011 ) for the density of occupied sites p{t) . Assuming the 



particles are spaced homogeneously in an infinite system 
we have 



dpjt) 
dt 



Pb{l-pd)p{tf{l-p{t))-Pdp{t). 



(2.2) 



The first term is the proliferation term and so is propor- 
tional to p^, the probability that the particle does not die 
[l—pd), the probability that the next randomly chosen site 
to give birth on is empty (1 — p) and finally the probability 
that it gives birth if this is the case, pj,. The second term 
represents particle annihilation and so is proportional to 
both p and pd, the probability that the chosen particle dies. 
Eq. (2.2) has three steady states, 



po = 0, p± = i(l±Jl — - 

2 \ V Pb{l-Pd) 



(2.3) 



^ We note that the birth rate is actually given by P(,(l — Pd) a-nd 
not Pi, only. 



For Apd > pb{l^Pd), P± are imaginary, resulting in po being 
the only real stationary state and so, here, extinction occurs 
in all circumstances. Keeping pb a constant from now on, 
we then have that our critical death rate is given by pd^ = 
Pb/{^ + Pb) which separates the active phase representing 
survival and the absorbing state of extinction. 

Clearly Eq. (2.2) is limited by the exclusion of diffusion 
and noise as well as the false assumption of a homogenous 
population density. We do however find at least good qual- 
itative support for our mean field analysis through numer- 
ical simulations. Fig. 1 a) shows the critical values of pd 
and Pb separating the regions with one and three station- 
ary states according to both the mean field equation and 
numerical simulations for 1, 2 and 3 dimensions. We see 
convincing agreement between our analytical and numer- 
ical results, particularly for higher dimensions. The MC 
simulations were carried out on an initially fully occupied 
lattice with linear sizes L — 1000, 32 and 10 for each di- 
mension respectively and we observed whether extinction 
occurred during 10^ time steps. For each birth rate, the 
simulation was repeated 500 times. If a single run survived, 
Pd was increased, whereas if extinction occurred in all runs, 
Pd was reduced. Using the same initial seed for the random 
number generator, an iterative procedure produced a crit- 
ical value with accuracy ±2^^^. This iterative procedure 
was then repeated 5 times with different seeds and the av- 
erage taken. Only a small number of repeats were needed 
since the largest variance of the values obtained was of the 
order of 10~^. From the figure we find that to 3 d.p. for 
Pb = 0.5, Pd, = 0.073 , 0.098 and 0.105 in 1,2 and 3 dimen- 
sions respectively. Due to the finite size of the lattices and 
the finite time used for the above simulations, the actual 
critical death rates are likely to differ slightly from those 
given and more acc urate techniques w ould have to be used 
to obtain them (see Hinrichsenl . 200d for examples of such 
techniques) . 

With pit = 0) = 1, as is increased, the steady-state 
population density decreases, becoming zero at pd, as 
shown in Fig. 1 b), marking the phase transition. We see 
that the steady state population density appears to change 
continuously in 1 dimension, whilst discontinuously in 2 
and 3 dimensions in agreement with the mean field results. 
If indeed this is the case, we call such phase transitions 
continuous and first-order respectively. In both cases, the 
phase transition is marked by a very rapid decrease in 
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Fig. 1. a) Parameter domain for the number of steady states with 
the points showing the parameters giving 2 steady states and b) the 
steady state population densities with pi, = 0.5 for the mean field 
(line), 1 (+), 2 (x) and 3 (•) dimensional simulations. 



population density. 

Briefly relating our model to biology, we note that the 
death rate for a given species may fluctuate for any number 
of reasons but it must certainly be true that at least the 
average value ofpd must be less than p^^ for the species to 
have ever been in existence. However, as we have shown, 
even a temporary increase in pd above pd^ , due to defor- 
estation or disease for example, will cause a very rapid and 
perhaps unrecoverable decrease in population. Extinctions 
however may also occur for reasons other than having a 
super-critical death rate. We investigate the roles of AUee 
effects and that of fluctuations in the next three sections 
where we examine simulations in the sub-critical or active 
phase and use the constant value pb = 0.5. 




Fig. 2. Phase diagram showing the critical values of separating the 

2 long-term outcomes of the system for different initial population 
density according to the mean field (line) and the 1 (-I-), 2 (x) and 

3 (•) dimensional MC simulations. 

3. Allee Effects 



One reason we observe a decline in population growth at 
low densities is due to individuals finding it harder to find 
a m ate. This is empirically known to occur in bot h plant 
(e.g. Aizen and Feinsingei . 1994} ) and animal (e.g. Lande . 
Il987l) populations. In our model, this aspect is incorporated 
by the fact that two individuals are required for reproduc- 
tion whereas an individual can die by itself. As density de- 
creases, each individual therefore finds it increasingly dif- 
ficult to find another for reproduction before they die. To 
examine this, we return to our mean field equation (2.2). 

It is easy to show that whereas p+ and po are stable 
stationary points of Eq. (2.2), p_ is unstable. Since in the 
active phase, po < p- < p+, any population whose density 
p{t) < p- will be driven to extinction by the dynamics of 
the system. In fact we find that for pd < Pd^ , 



Pit) 



for p{t) < p_ 
P+ for p{t) > p- 



as t 



(3.1) 



We test this numerically in 1, 2 and 3 spatial dimensions 
by finding the value oipd that separates the active and ab- 
sorbing states for different initial conditions. The MC sim- 
ulations were carried out and the critical death rate found 
iteratively in the same fashion as in Section 2. The results 
are shown in Fig. 2 and clearly show the importance of the 
initial population density for survival. The density depen- 
dence appears to increase with dimensionality, which we ex- 
pect, since two individuals meeting becomes progressively 
harder as the dimensionality of the system increases. 

The existence of this critical population density is highly 
significant to the conservation of species. It is clear that 
a sufficiently small population will not grow, regardless of 
how much space and resources are available. It also has 
repercussions if a population density were to suddenly de- 
crease due to disease or particularly harsh meteorological 
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Fig. 3. Log-log plot of the standard deviation of the population 
density versus the number of sites in the 1 (+), 2 (x) and 3 (•) 
dimensional systems. The hashed line has gradient -0.5 for the eye 
and indicates the power law behaviour. Insert: The fluctuations v.s. 
Pd dimensional case with the same symbol notation. 

conditions, for example. We examine this further in Section 
5 after examining the role of fluctuations. 




Fig. 4. How Ps varies with L for the 1 dimensional model with (from 
left to right), = 0.04, 0.05 and 0.06. Similar results are seen in 2 
and 3 dimensions. 

ing more susceptible to extinction. This is indeed observed 
in nature where long-lived species are known, in general, 
to have a hig her chance of survival than short-lived ones 
( Pimml . [l991h . 



4. Fluctuations 



5. A Decrease in Population Density 



We expect extinction due to fluctuations in the popu- 
lation density to occur when the order of the fluctuations 
approaches the mean population density. Empirically, de- 
mographic stochasticity (that is, chance events of mortality 
and rep roduction) is known to be greater in smaller popu- 
lations ( Lande et al than in larger ones. Population 
and habitat size are, on average, positively correlated and 
so, particularly with the existence of the critical population 
density, we expect extinction due to fluctuations to occur 
for smaller lattice sizes as has been s uggested by others (e.g. 
IPimml . I1991I : lEscudero et"all . l2004h . 

We see in Fig. 3 that, numerically, the fluctuations in the 
population density Ap decrease with the number of lattice 
sites N through a power law with exponent -0.50 in all di- 
mensions, which is what we would expect from the central 
limit theorem. Simulations were carried out for fixed = 
0.03 andpf, = 0.5 and the standard deviation obtained from 
5 X 10'^ surviving runs for each lattice size. The insert in 
Fig. 3 shows how the size of the fluctuations also increase as 
the critical point is approached. These larger fluctuations 
will also increase the probability of extinction as indicated 
in Fig. 4 where we examine the probability of survival Ps, 
that is, the probability that extinction has not occurred up 
to some time tm- We examine the 1 dimensional case only 
using three different values of pd with t„j = lO'^ and repeat 
the simulation 5 x 10* times for each lattice size. The fig- 
ures clearly show how the probability of survival increases 
with L, yet decreases as pd increases. Indeed, as pd^ is ap- 
proached, population density decreases and fluctuation size 
increases resulting in species with higher death rates be- 



Apart from the initial conditions, it is certainly conceiv- 
able that the population density could fall below the crit- 
ical value due to a reduction in population size. From Eq. 
(3.1) we expect that the population will survive only as 
long as p{t) > pc- We simulate this by increasing to 1 at 
some t = tk and then returning pd to what it was before, 
once a density ps has been reached. We examine this here 
in 2 dimensions with now reflective rather than the previ- 
ously used periodic boundary conditions. Qualitatively, all 
previous results have been very similar when using reflec- 
tive boundary conditions but here we want to increase this 
degree of realism in our model. 

For 2 dimensional simulations with an initial population 
density p{0) = 0.128, the critical death rate is 0.093 (3 d.p.) 
as shown in Fig. 2. So, for pd = 0.093, we would expect 
that if ps > 0.128, the population will survive, with the 
population density returning to what it was before, whereas 
for Ps < 0.128, extinction will occur. Due to the fluctuations 
in the population that occur in the simulations, we expect 
more of an increase in the likelihood of extinction as p — > 
p'^ rather than the definite survival/extinction result that 
the mean field predicts. Fig. 5 shows the results for ps = 
0.13, where we see that for those runs that did survive, the 
population density does indeed return to what it once was. 
We also see, as expected, that most of the runs did result 
in extinction. In fact the survival rate was 0.004. 

From Fig. 5 b) we observe that there is a time delay of 
approximately 40 time steps between the sudden decrease 
in population and when the survival probability begins to 
fall. Assuming a particle that dies the nth time it is picked, 
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Fig. 5. a) The average population density of the surviving runs only, 
b) The average population density of all the runs (solid line) and 
the survival probability P(t) (hashed line), i.e. the probability that 
extinction has not occurred up to time t. 
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Fig. 6. Plot showing the recovery of the population n{t) for the 
surviving runs only after a disease breakout at t = 3000 due to the 
re-sizing of the lattice. The lattice is returned to how it was originally 
at t = 6000 and the population recovers its original value. 
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Fig. 7. How the probability of survival changes with different reduc- 
tions in L, starting from L = 32. 



survives n— 1 time steps, it is easy to show that the expected 
Hfetime (in time steps) of an individual is given by (1 — 
Pd)/Pd- We therefore have a time delay of approximately 
four lifetimes (recall = 0.093) which, for a lot of species, 
is ample time to act. 

In order to prevent extinction in such a case, the popula- 
tion density must be increased beyond pc- This has impor- 
tant ecological implications since it shows that the proba- 
bility of extinction can be decreased, not only by increasing 
the population (which is of course not always possible), but 
also by a decrease in habitat size. 

To see whether this hypothesis holds, we simulate this 
again using = 0.093 but this time Ps = Pc = 0.128 so 
that the chance of survival is negligible. This time however, 
once the population density has been reduced, the area 
covered by the lattice is reduced by half. The organisms 
in the half that remains are left where they are, whereas 
those in the half that is removed are randomly placed in the 



remaining half. This then doubles the population density, 
bringing the population out of the sub-critical population 
density. Once the population has recovered and stabilised, 
the lattice size is returned to how it once was. The results 
are shown in Fig. 6 and clearly show the recovery of the 
population once the lattice size has been reduced. In fact, 
out of 1000 runs, the probability of survival rose from 0.003 
to 0.281. 

We expect there to be an optimal habitat reduction size 
- too large a reduction and the population will be in dan- 
ger from large fluctuations associated with smaller habitat 
sizes whereas too small a reduction and the density will 
not be increased sufficiently. We therefore plot in Fig. 7 the 
probability that the system does not go extinct up to some 
t = tmaxj Psi against the reduction in L, AL. We again use 
Pd = 0.093 and ps = 0.128. For small AL, changes little 
due to the density not being reduced enough, yet for larger 
AL, the larger fluctuations resulting from the smaller value 



of L also cause Pg to be small. 

Whilst reflective boundary conditions were used here, 
very similar results were obtained using periodic boundary 
conditions. In fact, with periodic boundary conditions, the 
probability of survival increased more significantly by the 
decrease in L due to the population being able to grow 
in two directions rather than in just one after the habitat 
size has been returned to what it once was. This of course 
could be achieved in reality by reducing the habitat from 
more than one direction. 

This model was proposed to represent how the area in 
which a population is found could be reduced in real-life. 
The species could be driven towards one end of the habi- 
tat with a boundary placed to prevent them leaving the 
desired area. This boundary could then be removed once 
the population has recovered. Clearly this is easier for 
larger, land-based animals but in principle, at least, could 
be achieved for all species. 
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Fig. 8. Snapshot of the output from a 2 dimensional lattice with 
L = 100. A value of p^j = 0.1 was used and the picture was taken at 
t = 600 when p = 0.2223. 



fact, for populations that are found in patches, the value 
of the density will depend very much on the scales used. 
The same is true of the MC results as shown in Figure 
8, where we see clear examples of clustering. In nature, 
species will cluster to varying degrees and hence the value 
Allee effects are cer t ainly observed in nature (iStephens and Swterfeadji ica.l population density will also vary and would 



6. Conclusions 



19991 : IPedersen et all . lioOlt ICvllenberg et al 



19971. et al) need to be estimated in each case. 



and have been studied with respect to extinction. Using a 
lattice model, we have observed Allee effects together with 
the role of fluctuations, with the advantage of being able 
to examine the effects of habitat size. Being able to model 
the population as a group of individuals which move, breed 
and die, rather than as a variable in an equation, has 
enabled us to gain a more realistic insight into how real 
populations behave. 

Rather than the clear-cut conclusions that deterministic 
models produce, conservationists often examine the proba- 
bility that a population will maintain itself without signif- 
icant demographic or genetic man ipulation for the foresee- 
able ecological future ( Soulel . 1987 ). In this spirit, for a suf- 
flciently large population density we have shown that the 
probability of survival does increase with habitat size due 
to the smaller fluctuations. However, far more important 
are the death rate and population density since if these fall 
on the wrong side of their critical values, extinction is al- 
most a certainty. 

Our findings are certainly significant for the design of 
habitats. The notion of a critical habitat size, mentioned in 
the Introduction, is misleading, since, it is certainly not true 
that for a fixed population size, the larger the habitat size 
the better. Regardless of the amount of space and resources 
available, a population will only grow if the density is above 
its critical value. We also proposed, in the last section, a 
method for greatly reducing the probability of extinction 
by reducing the habitat size once a species has become rare. 

Our notion of density has been that of the number of 
individuals per unit area. While we assumed this to be 
constant in space when deriving our mean field equation 
(2.2), clearly this will vary amongst real populations. In 



Compared to other stochastic models, we claim the use of 
lattice models gives a more realistic insight into the way 
in which real populations behave. We do still however, 
recognise the inaccuracies in our model and the difficulties 
in implementing the observations. We believe the model to 
be valid, to a greater or lesser degree to all species which 
rely on others for growth, perhaps particularly those who 
live alone yet sexually reproduce. In fact, due to the great 
variety of species, we have presented the above as ideas 
which may be of qualitative, rather than quantitative, 
relevance to conservation management. 
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